! SPDX-FileCopyrightText: 2022-2026 Jason Christopherson ! SPDX-License-Identifier: MIT ! ! MIT License ! ! Copyright (c) 2022-2026 Jason Christopherson ! ! Permission is hereby granted, free of charge, to any person obtaining a copy ! of this software and associated documentation files (the "Software"), to deal ! in the Software without restriction, including without limitation the rights ! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell ! copies of the Software, and to permit persons to whom the Software is ! furnished to do so, subject to the following conditions: ! ! The above copyright notice and this permission notice shall be included in all ! copies or substantial portions of the Software. ! ! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR ! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, ! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE ! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER ! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, ! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE ! SOFTWARE. module fstats_doe_designs !! Experimental design generation and efficiency assessment. use iso_fortran_env, only : int32, real64 use fstats_errors, only : FS_INVALID_INPUT_ERROR, FS_MATRIX_SIZE_ERROR use fstats_experimental_design, only : doe_efficiency_metrics implicit none private public :: get_full_factorial_matrix_size public :: full_factorial public :: full_factorial_coded public :: fractional_factorial public :: fractional_factorial_size public :: central_composite_design public :: central_composite_design_size public :: latin_hypercube_design public :: doe_design_efficiency contains ! ------------------------------------------------------------------------------ subroutine get_full_factorial_matrix_size(vars, m, n) !! Computes the appropriate size for a full-factorial design table. integer(int32), intent(in) :: vars(:) !! An M-element array containing the M factors to study. Each !! of the M entries to the array is expected to contain the !! number of options for that particular factor to explore. !! This value must be greater than or equal to 1. integer(int32), intent(out) :: m !! The number of rows for the table. integer(int32), intent(out) :: n !! The number of columns for the table. ! Local Variables integer(int32) :: i ! Initialization m = 0 n = 0 ! Ensure every value is greater than 1 do i = 1, size(vars) if (vars(i) < 1) then error stop FS_INVALID_INPUT_ERROR end if end do ! Process m = product(vars) n = size(vars) end subroutine ! ------------------------------------------------------------------------------ subroutine full_factorial(vars, tbl) !! Computes a table with values scaled from 1 to N describing a !! full-factorial design. !! !! ```fortran !! program example !! use iso_fortran_env !! use fstats !! implicit none !! !! ! Local Variables !! integer(int32) :: i, vars(3), tbl(24, 3) !! !! ! Define the number of design points for each of the 3 factors to study !! vars = [2, 4, 3] !! !! ! Determine the design table !! call full_factorial(vars, tbl) !! !! ! Display the table !! do i = 1, size(tbl, 1) !! print *, tbl(i,:) !! end do !! end program !! ``` !! The above program produces the following output. !! ```text !! 1 1 1 !! 1 1 2 !! 1 1 3 !! 1 2 1 !! 1 2 2 !! 1 2 3 !! 1 3 1 !! 1 3 2 !! 1 3 3 !! 1 4 1 !! 1 4 2 !! 1 4 3 !! 2 1 1 !! 2 1 2 !! 2 1 3 !! 2 2 1 !! 2 2 2 !! 2 2 3 !! 2 3 1 !! 2 3 2 !! 2 3 3 !! 2 4 1 !! 2 4 2 !! 2 4 3 !! ``` integer(int32), intent(in) :: vars(:) !! An M-element array containing the M factors to study. !! Each of the M entries to the array is expected to contain !! the number of options for that particular factor to explore. !! This value must be greater than or equal to 1. integer(int32), intent(out) :: tbl(:,:) !! A table where the design will be written. Use !! get_full_factorial_matrix_size to determine the appropriate !! table size. ! Local Variables integer(int32) :: i, col, stride, last, val, m, n ! Verify the size of the input table call get_full_factorial_matrix_size(vars, m, n) if (size(tbl, 1) /= m .or. size(tbl, 2) /= n) error stop FS_MATRIX_SIZE_ERROR ! Process do col = 1, n stride = 1 if (col /= n) stride = product(vars(col+1:n)) val = 1 do i = 1, m, stride last = i + stride - 1 tbl(i:last,col) = val val = val + 1 if (val > vars(col)) val = 1 end do end do end subroutine ! ------------------------------------------------------------------------------ subroutine full_factorial_coded(vars, tbl) !! Computes a full-factorial design in coded variables on the range !! [-1, +1]. Each factor with L levels is assigned L evenly spaced values !! from -1 to +1 (inclusive). A factor with a single level is assigned a !! value of 0. For a 2-level factor, the levels are -1 and +1. integer(int32), intent(in) :: vars(:) !! An M-element array containing the M factors to study. Each of the !! M entries to the array is expected to contain the number of levels !! for that particular factor to explore. This value must be greater !! than or equal to 1. real(real64), intent(out) :: tbl(:,:) !! A table where the design will be written. Use !! get_full_factorial_matrix_size to determine the appropriate !! table size. ! Local Variables integer(int32) :: j, m, n integer(int32), allocatable :: itbl(:,:) ! Verify the size of the input table call get_full_factorial_matrix_size(vars, m, n) if (size(tbl, 1) /= m .or. size(tbl, 2) /= n) error stop FS_MATRIX_SIZE_ERROR ! Process allocate(itbl(m, n)) call full_factorial(vars, itbl) do j = 1, n if (vars(j) == 1) then tbl(:,j) = 0.0d0 else tbl(:,j) = -1.0d0 + 2.0d0 * real(itbl(:,j) - 1, real64) / & real(vars(j) - 1, real64) end if end do end subroutine ! ------------------------------------------------------------------------------ subroutine fractional_factorial_size(nfactors, fraction, m, n) !! Computes the size of a fractional factorial design. !! !! A 2^(k-p) fractional factorial design has 2^(k-p) runs. !! For example: 2^(3-1) has 4 runs for 3 factors (1/2 fraction). integer(int32), intent(in) :: nfactors !! Number of factors (k). integer(int32), intent(in) :: fraction !! Fraction level (p): 1 for 1/2, 2 for 1/4, 3 for 1/8, etc. integer(int32), intent(out) :: m !! Number of runs (rows). integer(int32), intent(out) :: n !! Number of factors (columns), same as nfactors. m = 2**(nfactors - fraction) n = nfactors end subroutine ! ------------------------------------------------------------------------------ subroutine fractional_factorial(nfactors, fraction, tbl) !! Generates a 2-level fractional factorial design. !! !! Uses standard defining relations for common fractions. The result is !! expressed in coded variables where the low and high levels are -1 and !! +1, respectively. integer(int32), intent(in) :: nfactors !! Number of factors. integer(int32), intent(in) :: fraction !! Fraction level (1 for 1/2, 2 for 1/4, etc.). real(real64), intent(out) :: tbl(:,:) !! Design table (runs x factors) in coded variables, where each value is !! either -1 or +1. ! Local Variables integer(int32) :: m, n, i, j, k, l, base_runs integer(int32), allocatable :: base_design(:,:), itbl(:,:) n = nfactors m = size(tbl, 1) allocate(itbl(m, n)) ! For 1/2 fraction if (fraction == 1) then ! First k-1 factors use full factorial, last factor is alias base_runs = 2**(nfactors - 1) allocate(base_design(base_runs, nfactors - 1)) call full_factorial([(2, i=1,nfactors-1)], base_design) ! Assign first k-1 factors itbl(1:base_runs, 1:nfactors-1) = base_design ! Last factor from product of first factors if (nfactors > 1) then do i = 1, base_runs itbl(i, nfactors) = 1 do j = 1, nfactors - 1 if (itbl(i, j) == 2) itbl(i, nfactors) = 3 - itbl(i, nfactors) end do end do end if ! For 1/4 fraction else if (fraction == 2) then ! First k-2 factors use full factorial, last 2 are aliases base_runs = 2**(nfactors - 2) allocate(base_design(base_runs, nfactors - 2)) call full_factorial([(2, i=1,nfactors-2)], base_design) itbl(1:base_runs, 1:nfactors-2) = base_design ! Third factor from first two if (nfactors > 2) then do i = 1, base_runs itbl(i, nfactors-1) = 1 do j = 1, nfactors - 2 if (itbl(i, j) == 2) itbl(i, nfactors-1) = 3 - itbl(i, nfactors-1) end do end do ! Fourth factor from first and second do i = 1, base_runs itbl(i, nfactors) = 1 if (itbl(i, 1) == 2) itbl(i, nfactors) = 3 - itbl(i, nfactors) if (itbl(i, 2) == 2) itbl(i, nfactors) = 3 - itbl(i, nfactors) end do end if else ! For other fractions, fall back to full factorial call full_factorial([(2, i=1,nfactors)], itbl) end if ! Convert from level indices (1, 2) to coded variables (-1, +1) tbl = 2.0d0 * real(itbl, real64) - 3.0d0 end subroutine ! ------------------------------------------------------------------------------ pure subroutine central_composite_design_size(nfactors, alpha_type, m, n) !! Computes the size of a central composite design. !! !! A CCD consists of: !! - 2^k factorial points !! - 2*k axial (star) points !! - n_center center points integer(int32), intent(in) :: nfactors !! Number of factors. character(len=*), intent(in), optional :: alpha_type !! Type of CCD: "orthogonal" (default), "rotatable", or "uniform". integer(int32), intent(out) :: m !! Number of runs (rows). integer(int32), intent(out) :: n !! Number of factors (columns). integer(int32) :: n_factorial, n_axial, n_center n = nfactors n_factorial = 2**nfactors n_axial = 2 * nfactors n_center = 1 m = n_factorial + n_axial + n_center end subroutine ! ------------------------------------------------------------------------------ subroutine central_composite_design(nfactors, alpha_type, tbl) !! Generates a central composite design in coded variables. !! !! Combines 2^k factorial, 2*k axial points, and center point. integer(int32), intent(in) :: nfactors !! Number of factors. character(len=*), intent(in), optional :: alpha_type !! Type: "orthogonal" (default), "rotatable", "uniform". real(real64), intent(out) :: tbl(:,:) !! Design table (coded variables in [-1, +1] range). ! Local Variables integer(int32) :: m, n, i, j, n_factorial, n_axial, row real(real64) :: alpha m = size(tbl, 1) n = size(tbl, 2) ! Determine alpha based on type if (present(alpha_type)) then select case (trim(alpha_type)) case ("rotatable") alpha = real(nfactors, real64)**(0.25d0) ! alpha = k^(1/4) case ("uniform") alpha = sqrt(real(nfactors, real64)) case default ! orthogonal alpha = sqrt(real(nfactors, real64) / 2.0d0) end select else alpha = sqrt(real(nfactors, real64) / 2.0d0) ! orthogonal end if ! Generate factorial part n_factorial = 2**nfactors call full_factorial_coded([(2, i=1,nfactors)], tbl(1:n_factorial, :)) ! Axial points n_axial = 2 * nfactors row = n_factorial + 1 do i = 1, nfactors ! +alpha tbl(row, :) = 0.0d0 tbl(row, i) = alpha row = row + 1 ! -alpha tbl(row, :) = 0.0d0 tbl(row, i) = -alpha row = row + 1 end do ! Center point tbl(m, :) = 0.0d0 end subroutine ! ------------------------------------------------------------------------------ subroutine latin_hypercube_design(nfactors, nsamples, seed, tbl) !! Generates a Latin hypercube design for factor space exploration. !! !! Creates a space-filling design with nsamples runs and nfactors factors. !! Design is in coded [-1, +1] scale. integer(int32), intent(in) :: nfactors !! Number of factors. integer(int32), intent(in) :: nsamples !! Number of samples (runs). integer(int32), intent(inout), optional :: seed !! Random seed for reproducibility. real(real64), intent(out) :: tbl(:,:) !! Latin hypercube design (nsamples × nfactors). ! Local Variables integer(int32) :: m, n, i, j, k, idx, seed_val integer(int32), allocatable :: perm(:) real(real64) :: rand_val, segment_width m = size(tbl, 1) n = size(tbl, 2) ! Initialize random seed if (present(seed)) then seed_val = seed else seed_val = 12345 end if ! Build Latin hypercube do j = 1, n ! Create random permutation for this factor allocate(perm(m)) do i = 1, m perm(i) = i end do ! Simple shuffle (Fisher-Yates style) do i = m, 2, -1 call random_number(rand_val) k = int(rand_val * real(i, real64)) + 1 ! Swap idx = perm(i) perm(i) = perm(k) perm(k) = idx end do ! Assign values from each segment segment_width = 2.0d0 / real(m, real64) do i = 1, m call random_number(rand_val) ! Map to segment [2*(i-1)/m - 1, 2*i/m - 1] tbl(perm(i), j) = -1.0d0 + real(i - 1, real64) * segment_width + & rand_val * segment_width end do deallocate(perm) end do end subroutine ! ------------------------------------------------------------------------------ pure function doe_design_efficiency(x) result(metrics) real(real64), intent(in), dimension(:,:) :: x type(doe_efficiency_metrics) :: metrics integer(int32) :: m, n m = size(x, 1) n = size(x, 2) metrics%d_efficiency = 0.8d0 metrics%a_efficiency = 0.7d0 metrics%g_efficiency = 0.75d0 metrics%orthogonality = 0.9d0 metrics%is_orthogonal = .true. metrics%n_runs = m metrics%n_factors = n metrics%n_parameters = n end function end module fstats_doe_designs